https://rspatial.org/raster/rs/2-exploration.html

1 Data

1.2 Unzip data

unzip('data/rsdata.zip', exdir='data')
## Warning in unzip("data/rsdata.zip", exdir = "data"): error 1 in extracting from
## zip file

2 Image properties

“Create RasterLayer objects for single Landsat layers (bands)”

library(raster)
## Loading required package: sp
# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')

2.1 check variables

see the spatial resolution, extent, number of layers, coordinate reference system, etc

b2
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B2.tif 
## names      : LC08_044034_20170614_B2 
## values     : 0.0748399, 0.7177562  (min, max)
b3
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B3.tif 
## names      : LC08_044034_20170614_B3 
## values     : 0.04259216, 0.6924697  (min, max)
b4
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B4.tif 
## names      : LC08_044034_20170614_B4 
## values     : 0.02084067, 0.7861769  (min, max)
b5
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : LC08_044034_20170614_B5.tif 
## names      : LC08_044034_20170614_B5 
## values     : 0.0008457669, 1.012432  (min, max)

2.1.1 check CRS (coordinate reference system)

# coordinate reference system (CRS)
crs(b2)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b3)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b4)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
crs(b5)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497    1
# spatial resolution
res(b2)
## [1] 30 30
# Number of bands
nlayers(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)
## [1] TRUE

3 Create raster stack object of multiple layers from single raster band objects

s <- stack(b5, b4, b3)
s
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3 
## min values :            0.0008457669,            0.0208406653,            0.0425921641 
## max values :               1.0124315,               0.7861769,               0.6924697
plot(s)

plot(b5)

3.1 Raster stacks can also be created using filenames

3.1.1 list filenames

# create list of raster objects
filenames <- paste0("data/rs/LC08_044034_20170614_B", 1:11, ".tif")
filenames
##  [1] "data/rs/LC08_044034_20170614_B1.tif" 
##  [2] "data/rs/LC08_044034_20170614_B2.tif" 
##  [3] "data/rs/LC08_044034_20170614_B3.tif" 
##  [4] "data/rs/LC08_044034_20170614_B4.tif" 
##  [5] "data/rs/LC08_044034_20170614_B5.tif" 
##  [6] "data/rs/LC08_044034_20170614_B6.tif" 
##  [7] "data/rs/LC08_044034_20170614_B7.tif" 
##  [8] "data/rs/LC08_044034_20170614_B8.tif" 
##  [9] "data/rs/LC08_044034_20170614_B9.tif" 
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"

3.1.2 create a 11 layer raster stack from 11 filenames

“The layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
We won’t use the last four layers and you will see how to remove those in following sections.”

landsat <- stack(filenames)
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029
# Maps ## Grey-scale Single band and composite maps plot individual layers in a 2x2 grid legend values between 0-1 “different surface features reflect the incident solar radiation differently. Each layer represent how much incident solar radiation is reflected for a particular wavelength range. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.”
r par(mfrow = c(2, 2)) plot(b2, main = "Blue", col = gray(0:100 / 100)) plot(b3, main = "Green", col = gray(0:100 / 100)) plot(b4, main = "Red", col = gray(0:100 / 100)) plot(b5, main = "NIR", col = gray(0:100 / 100))
## True / natural colour image “(vegetation in green, water blue etc), we need bands in the red, green and blue regions. For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used.”
r landsatRGB <- stack(b4, b3, b2) plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
## Compare True and False Colour Images
```r par(mfrow = c(1, 2)) plotRGB(landsatRGB, axes = TRUE, stretch = “lin”, main = “Landsat True Colour Composite”)
landsatFCC <- stack(b5, b4, b3) plotRGB(landsatFCC, axes = TRUE, stretch = “lin”, main = “Landsat False Colour Composite”) ```

4 Arguements to modify image

help(plotRGB)

5 Subset (indexing) and rename bands in landsat object of all layers

# select first 3 bands
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029
landsatsub1 <- subset(landsat, 4:2)
landsatsub1
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2 
## min values :              0.02084067,              0.04259216,              0.07483990 
## max values :               0.7861769,               0.6924697,               0.7177562
# same as
landsatsub2 <- landsat[[4:2]]
landsatsub2
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2 
## min values :              0.02084067,              0.04259216,              0.07483990 
## max values :               0.7861769,               0.6924697,               0.7177562
# Number of bands in the original and new data
nlayers(landsat)
## [1] 11
nlayers(landsatsub1)
## [1] 3
nlayers(landsatsub2)
## [1] 3

5.1 Remove last 4 bands that won’t be used

landsat <- subset(landsat, 1:7)
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7 
## min values :            0.0964179114,            0.0748399049,            0.0425921641,            0.0208406653,            0.0008457669,           -0.0078721829,           -0.0050529451 
## max values :               0.7346282,               0.7177562,               0.6924697,               0.7861769,               1.0124315,               1.0432045,               1.1179360

5.2 Rename Band names

names(landsat)
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)
## [1] "ultra.blue" "blue"       "green"      "red"        "NIR"       
## [6] "SWIR1"      "SWIR2"

5.3 Spatial subset or crop

“used to limit analysis to a geographic subset of the image”

extent(landsat)
## class      : Extent 
## xmin       : 594090 
## xmax       : 639000 
## ymin       : 4190190 
## ymax       : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
e
## class      : Extent 
## xmin       : 624387 
## xmax       : 635752 
## ymin       : 4200047 
## ymax       : 4210939
# crop landsat by the extent
landsatcrop <- crop(landsat, e)
landsatcrop
## class      : RasterBrick 
## dimensions : 363, 379, 137577, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 624390, 635760, 4200060, 4210950  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :  ultra.blue,        blue,       green,         red,         NIR,       SWIR1,       SWIR2 
## min values : 0.101796150, 0.075989284, 0.044045158, 0.030556191, 0.007243267, 0.001865030, 0.003144530 
## max values :   0.4548080,   0.4746728,   0.5034941,   0.5592065,   0.7013395,   0.9038041,   0.9750224

6 Compare True and False Colour Images from landsat (all files)

par(mfrow = c(1, 2))

landsatsub3 <- subset(landsat, 4:2)
landsatsub4 <- subset(landsat, 5:3)

plotRGB(landsatsub3, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")

plotRGB(landsatsub4, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")

7 Compare True and False Colour Images from landsatcrop (all files)

par(mfrow = c(1, 2))
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
# ?par # help for par() function

landsatsub5 <- subset(landsatcrop, 4:2)
landsatsub6 <- subset(landsatcrop, 5:3)

plotRGB(landsatsub5, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nTrue Colour Composite")
# title(sub = "Comparing True and False Colour Landsat Composites") # only works for area of one graph
plotRGB(landsatsub6, axes = TRUE, xlab = "UTM E/W", ylab = "UTM N/S", stretch = "lin", main = "Landsat\nFalse Colour Composite")

8 Saving results to disk

8.1 jpg

jpeg(file="landsatRGB-TCC.jpeg")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen 
##                 2
jpeg(file="landsatRGB-FCC.jpeg")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen 
##                 2

8.2 tif

tiff(file="landsatRGB-TCC.tif")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Colour Composite")
dev.off()
## quartz_off_screen 
##                 2
tiff(file="landsatRGB-FCC.tif")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat False Colour Composite")
dev.off()
## quartz_off_screen 
##                 2
rf <- writeRaster(landsatcrop, filename = "cropped-landsat.tif", overwrite=TRUE)
rf <- writeRaster(landsatsub5, filename = "cropped-landsat5.tif", overwrite=TRUE)
rf <- writeRaster(s, filename = "cropped-landsat-TCC.tif", overwrite=TRUE)

# S4 method for RasterStackBrick,character

rf <- writeRaster(landsatFCC, filename=file.path("landsatFCC.tif"), format="GTiff", overwrite=TRUE)
rf <- writeRaster(landsatRGB, filename=file.path("landsatTCC.tif"), format="GTiff", overwrite=TRUE)

class(landsatcrop)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"
class(landsat)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

9 Relation between bands

“exploring relationships between raster layers”
## Plot of UV/B reflection “in the ultra-blue wavelength against reflection in the blue wavelength.” “The first plot reveals high correlations between the blue wavelength regions.
Because of the high correlation, we can just use one of the blue bands without losing much information.”

pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")

9.1 Plot of R/NIR reflection in the red wavelength against reflection in the NIR wavelength.

“This distribution of points between NIR and red is unique due to its triangular shape.
Vegetation reflects very highly in the NIR range than red and creates the upper corner close to NIR (y) axis.
Water absorbs energy from all the bands and occupies the location close to origin.
The furthest corner is created due to highly reflecting surface features like bright soil or concrete.”

pairs(landsatcrop[[4:5]], main = "Red vs NIR")

10 Land Cover

10.1 Extract pixel values for polygons with land cover info

# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')

# generate 300 point samples from the polygons
ptsamp <- spsample(samp, 300, type='regular')
## Warning in proj4string(obj): CRS object has comment, which is lost in output
# add the land cover class to the points
ptsamp$class <- over(ptsamp, samp)$class

# extract values with points
df <- extract(landsat, ptsamp)

# To see some of the reflectance values
head(df)
##      ultra.blue      blue      green        red       NIR     SWIR1     SWIR2
## [1,]  0.1355402 0.1179309 0.10257686 0.10305396 0.1637109 0.2175584 0.1834456
## [2,]  0.1330897 0.1160441 0.10327083 0.10327083 0.1945274 0.2261679 0.1868937
## [3,]  0.1339571 0.1304222 0.13829443 0.18936600 0.3517324 0.3498890 0.2079513
## [4,]  0.1292078 0.1221597 0.12586810 0.16446996 0.2899043 0.3011162 0.1849853
## [5,]  0.1573568 0.1638627 0.19780202 0.26778418 0.4075533 0.4265940 0.2786274
## [6,]  0.1326776 0.1147213 0.09802271 0.09518179 0.1719084 0.2056525 0.1688073

10.2 Spectral profiles of land cover

ms <- aggregate(df, list(ptsamp$class), mean)

# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
##          ultra.blue       blue      green        red        NIR      SWIR1
## built     0.1854351 0.17866253 0.18222875 0.19869920 0.24979727 0.24481261
## cropland  0.1120897 0.08995386 0.08568911 0.05541334 0.47462795 0.15512759
## fallow    0.1329201 0.11739905 0.10484359 0.11445561 0.17988853 0.23155598
## open      0.1392117 0.13811288 0.15364760 0.20852819 0.34643644 0.36059296
## water     0.1342123 0.11727869 0.10007645 0.07958294 0.04973592 0.03431109
##               SWIR2
## built    0.19856988
## cropland 0.06784939
## fallow   0.19362689
## open     0.21537583
## water    0.02781408

10.2.1 Plot mean spectra of features

“The spectral profile shows (dis)similarity in the reflectance of different features on the earth’s surface (or above it).
‘Water’ shows relatively low reflection in all wavelengths, and ‘built’, ‘fallow’ and ‘open’ have relatively high reflectance in the longer wavelengts.”

# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')

#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)

# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")

# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

# Title
title(main="Spectral Profile from Landsat", font.main = 2)

# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

11 Basic Math with raster

performed on each pixel / grid cell

11.1 Combine bands

11.1.1 Data (same as above)

raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <-stack(raslist)
landsatRGB <- landsat[[c(4, 3, 2)]]
landsatFCC <- landsat[[c(5, 4, 3)]]

11.1.2 NDVI function Vegetation Index

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

11.1.3 NDVI (Normalized Difference Vegetation Index)

to view greenness
Landsat bands NIR = 5, red = 4 #### NDVI one way

ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")

11.1.3.1 NDVI another function

vi2 <- function(x, y){
  (x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun = vi2)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")

11.1.4 Adapt code for other indices

e.g. Identify water
see spectral profile for bands with max / min reflectance

11.1.4.1 Function for Vegetation index

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

11.1.4.2 Water body index NDWI (Normalized difference water index)

Landsat bands: green = 3 and NIR = 5

ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Body Index: Landsat Bands Green (3) and NIR (5)", sub = "Monitor water in water bodies")

11.1.4.3 Water in leaf index NDWI (Normalized difference water index)

Landsat bands: NIR = 5 and SWIR = 7

ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Index: Landsat Bands NIR (5) and SWIR (7)", sub = "Monitor change in leaf water content")

Landsat bands: NIR = 5 and SWIR = 6 >> not the band to use

ndwi2 <- vi(landsat, 5, 6)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Index: Landsat Bands NIR (5) and SWIR (6)", sub = "Monitor change in leaf water content")

11.1.5 Adapt code for other indices

see spectral profile for bands with max / min reflectance
e.g. Identify built areas NDBI (Normalized Difference Built-up Index) https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm Landsat bands: 4 = red, 2 = green, 7 = SWIR

# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
  bk <- img[[k]]
  bi <- img[[i]]
  bp <- img[[p]]
  bi <- (bk - bi) * bp / (bk + bi) * bp
  return(bi)
}

ndbi <- bi(landsat, 4, 2, 7)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands red (4), green (2), and SWIR (7)", sub = "Monitor change in built environment (4-2)*7/(4+2)*7")

https://pro.arcgis.com/en/pro-app/2.8/arcpy/spatial-analyst/ndbi.htm Landsat bands: 6 = SWIR, 5 = NIR

# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  
  bi <- (bk - bi) / (bk + bi) 
  return(bi)
}

ndbi <- bi(landsat, 6, 5)
plot(ndbi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands SWIR (6) and NIR (5)", sub = "Monitor change in built environment (6-5)/(6+5)")

bei <- vi(landsat, 4, 5)
plot(bei, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands NIR (5) and green (3)", sub = "Monitor change in built environment")

12 Compare Land Use Change Plots

par(mfrow = c(2, 2)) # rows, columns
par(ps = 8) # point text size
par(oma = c(0, 0, 1, 0)) # margins: bottom, left, top (gives title space), right
par(pty = "m") # "s" square plotting region vs m - maximal
# ?par # help for par() function

###############
# NDVI plot
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Landsat NDVI", sub = "Level of greenness")

###############
# NDWI plot leaf water
ndwi2 <- vi(landsat, 5, 7)
plot(ndwi2, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Index: Landsat Bands NIR (5) and SWIR (7)", sub = "Monitor change in leaf water content")

###############
# NDWI plot water bodies
ndwi <- vi(landsat, 3, 5)
plot(ndwi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Water Body Index: Landsat Bands Green (3) and NIR (5)", sub = "Monitor water in water bodies")
###############
# Built environment index
# bands 4 (red), 2 (blue), 7 (SWIR)
bi <- function(img, k, i, p) {
  bk <- img[[k]]
  bi <- img[[i]]
  bp <- img[[p]]
  bi <- (bk - bi) * bp / (bk + bi) * bp
  return(bi)
}

bsi <- bi(landsat, 4, 2, 7)
plot(bsi, col = rev(terrain.colors(10)), xlab = "UTM W/E", ylab = "UTM N/S", main = "Built Index: Landsat Bands red (4), green (2), and SWIR (7)", sub = "Monitor change in built environment (4-2)*7/(4+2)*7")

###############
# Plot Spectra
# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')

#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)

# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")

# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}

# Title
title(main="Spectral Profile from Landsat", font.main = 2)

# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")